#ifndef MATRIX_OP_H_INCLUDED
#define MATRIX_OP_H_INCLUDED

typedef struct eigensystem
{
    float **eigenvalues;
    float **eigenvector;
}eigen;

float **transpose(float **mat, int h, int w)
{
    float **res=NULL;
    int row, col;
    res=alloc2Dfloat(w,h);
	
    for(row=0;row<h;row++)
    {
        for(col=0;col<w; col++)
        {
           
		 res[col][row]=mat[row][col];
        }
    }
    return res;
}

float **add_mat(float **mat1, float **mat2, int h, int w)
{
    float **res;
    int row, col;
    res=alloc2Dfloat(h, w);
    for(row=0;row<h; row++)
    {
        for(col=0;col<w;col++)
        {
            res[row][col]=mat1[row][col]+mat2[row][col];
        }
    }
    return res;
}

float **multiply_mat(float **mat1, float **mat2, int h1, int w1, int h2, int w2)
{
    float **res, temp;
    int row, col, r_ind;
    res=alloc2Dfloat(h1, w2);
    for(row=0;row<h1;row++)
    {
      for(col=0;col<w2;col++)
      {
          temp=0.0;
          for(r_ind=0;r_ind<w1;r_ind++)
          {
              temp=temp+mat1[row][r_ind]*mat2[r_ind][col];
          }
          res[row][col]=temp;
      }
    }
   return res;

}

float **inv(float **matA, int n)
{

    int **matI,i,j,maxi,k;
    float **matB,temp,temp1,**matres;
    matI=(int**)calloc(n,sizeof(int**));
    for(i=0;i<n;i++)
    {
               matI[i]=(int*)calloc(n,sizeof(int**));

    }
     for(i=0;i<n;i++)
    {
                    for(j=0;j<n;j++)
                    {
                        matI[i][j]=i==j ? 1 : 0;
                    }
    }
       matB=(float**)calloc(n,sizeof(float**));
    for(i=0;i<2*n;i++)
    {
               matB[i]=(float*)calloc(2*n,sizeof(float*));

    }
    for(i=0;i<n;i++)
    {
                    for(j=0;j<n;j++)
                    {
                        matB[i][j]=matA[i][j];
                        matB[i][j+n]=matI[i][j];
                    }

   }

  for(i=0;i<n;i++)
  {
                  maxi=i;
                  temp=matB[i][i];
                  for(j=i+1;j<n;j++)
                  {
                                    if(fabs(matB[j][i])>fabs(temp))
                                    {
                                                                   maxi=j;
                                                                   temp=matB[j][i];
                                    }
                  }
       	for(j=0;j<2*n;j++)
		{
			temp1=matB[i][j];
			matB[i][j]=matB[maxi][j];
			matB[maxi][j]=temp1;
		}
	temp=matB[i][i];
    for(j=0;j<2*n;j++)
		{
		  if(matB[i][j]>0 && temp>0)
		  {
			matB[i][j]=(float)matB[i][j]/temp;
		  }
		  else
		  {
		    matB[i][j]=0;
		  }

		}
  for(k=0;k<n;k++)
		{
			if(k==i)
				continue;
			temp=matB[k][i];
			for(j=0;j<2*n;j++)
			{
				matB[k][j]=matB[k][j]-(temp*matB[i][j]);


			}
		}
}
matres=(float**)calloc(n,sizeof(float**));
    for(i=0;i<n;i++)
    {
               matres[i]=(float*)calloc(n,sizeof(float**));

    }
 for(i=0;i<n;i++)
 {
                 for(j=0;j<n;j++)
                 {
                                 matres[i][j]=matB[i][j+n];
                 }
}

 return matres;

}
float **Identity_mat(int n)
{
    float **V;
    int row;
    V=alloc2Dfloat(n,n);
    for(row=0;row<n;row++)
    {
        V[row][row]=1.0;
    }
    return V;
}

float **sym_schur(float **A, int n, int p, int q)
{
         float **J, c, s,r, t;
         J=Identity_mat(n);
         if(A[p][q]!=0)
         {
             r=(A[q][q]-A[p][p])/(2*A[p][q]);
             if(r>0)
              {
                  t=1/(r+sqrt((float)(1+(r*r))));
              }
              else
              {
                  t=-1/(-r+sqrt((float)(1+(r*r))));
              }
              c=1/(sqrt(1+t*t));
              s=t*c;
         }
         else
         {
             c=1;
             s=0;
         }
       J[p][p]=c;
       J[p][q]=s;
       J[q][p]=-s;
       J[q][q]=c;
 return J;
}
void copy_mat(float **A, float **mat, int m, int n)
{
    int r, c;
    for(r=0;r<m;r++)
    {
        for(c=0;c<n;c++)
        {
            A[r][c]=mat[r][c];
        }
    }
    
}
/*eigen find_eigenValue_eigenvector(float **mat, int n)
{
    float **V, **J, **A, **J_trans, **temp;
    int p, q;
    eigen e;
    V=Identity_mat(n);
    A=alloc2Dfloat(n,n);
    A=copy(mat,A, n, n);


    for(p=0;p<n-1;p++)
    {
        for(q=p+1;q<n;q++)
        {
          J=sym_schur(A,n,p,q);
          J_trans=transpose(J,n,n);
          temp=multiply_mat(J_trans,A,n,n,n,n);
          A=multiply_mat(temp,J,n,n,n,n);
          V=multiply_mat(V, J, n,n,n,n);
        }
    }
    e.eigenvalues=A;
    e.eigenvector=V;
    return e;
}
float *get_eigenvectors(float **mat, int n)
{
    float *w, *temp;
    int *ind, r;
    eigen e;
    temp=alloc1Dfloat(n);

    e=find_eigenValue_eigenvector(mat, n);
    for (r=0;r<n;r++)
    {
        temp[r]=e.eigenvalues[r][r];
    }
    ind=sortedindx(temp, n);
    w=alloc1Dfloat(n);
    for(r=0;r<n;r++)
    {
        w[r]=e.eigenvector[ind[0]][r];
    }
    return w;
}

float **diagonalize_matrix(float *w, int n)
{
  float **mat_new;
  int i, j;
  mat_new=alloc2Dfloat(n,n);
  for(i=0;i<n;i++)
    {
      for(j=0;j<n;j++)
	{
	  if(i==j)
	    {
	      mat_new[i][j]=w[i];
	  
	    }
	  else
	    {
	      mat_new[i][j]=0;
	    }
        }
    }
  return mat_new;
}

float **augumen_mat_val(float **mat, int m,int n, int aug_val)
{
 
  int i;
  for(i=0;i<m;i++)
    {
      mat[i]=(float*)realloc(mat[i],n+1*sizeof(float)); 
      mat[i][n]=aug_val;
    }
  
  return mat;      

}*/

float mean2D(float **matA, int height, int width)
{
	float mean=0.0;
	int r_indx, c_indx;
	for(r_indx=0;r_indx<height; r_indx++)
	{
		for(c_indx=0;c_indx<width;c_indx++)
		{
			mean=mean+matA[r_indx][c_indx];
		}
		
	}
	mean=mean/(height*width);
	return mean;
	
}

float **vec2mat(float *vec, int m)
{
	
	float **mat;
	int row;
	mat=alloc2Dfloat(m, 1);
	for(row=0;row<m; row++)
	{
		mat[row][0]=vec[row];
	}
	return mat;
}




#else
#endif
